#include "data_structures.hpp"
#include <iostream>
#include <iomanip>
#include <numeric>

namespace EddyDetection {

nlohmann::json GeoJSONExporter::createFeature(const EddyContour& eddy,
                                               const std::string& method_name,
                                               double depth_m) {
    nlohmann::json feature;
    
    feature["type"] = "Feature";
    
    // Creating Geometry
    nlohmann::json geometry;
    geometry["type"] = "Polygon";
    
    // Create coordinate array
    nlohmann::json coordinates = nlohmann::json::array();
    nlohmann::json ring = nlohmann::json::array();
    
    // Add boundary points
    for (const auto& point : eddy.boundary_points) {
        nlohmann::json coord = nlohmann::json::array();
        coord.push_back(point[0]);
        coord.push_back(point[1]);
        ring.push_back(coord);
    }
    
    // Ensure polygon closure
    if (!eddy.boundary_points.empty()) {
        const auto& first_point = eddy.boundary_points[0];
        const auto& last_point = eddy.boundary_points.back();
        
        if (std::abs(first_point[0] - last_point[0]) > 1e-6 ||
            std::abs(first_point[1] - last_point[1]) > 1e-6) {
            nlohmann::json coord = nlohmann::json::array();
            coord.push_back(first_point[0]);
            coord.push_back(first_point[1]);
            ring.push_back(coord);
        }
    }
    
    coordinates.push_back(ring);
    geometry["coordinates"] = coordinates;
    feature["geometry"] = geometry;
    
    // Creating properties
    nlohmann::json properties;
    properties["eddy_id"] = eddy.id;
    properties["detection_method"] = method_name;
    properties["analysis_depth_m"] = std::round(depth_m * 100.0) / 100.0;
    properties["type"] = eddy.type;
    properties["thermal_type"] = eddy.thermal_type;
    properties["center_lat"] = std::round(eddy.center[0] * 10000.0) / 10000.0;
    properties["center_lon"] = std::round(eddy.center[1] * 10000.0) / 10000.0;
    properties["area_km2"] = std::round(eddy.area_km2 * 100.0) / 100.0;
    properties["mean_vorticity"] = std::round(eddy.mean_vorticity * 100000000.0) / 100000000.0;
    properties["temp_anomaly_C"] = std::round(eddy.temp_anomaly * 1000.0) / 1000.0;
    
    feature["properties"] = properties;
    
    return feature;
}

nlohmann::json GeoJSONExporter::createGeoJSONObject(const std::vector<EddyContour>& eddies,
                                                     const std::string& method_name,
                                                     double depth_m) {
    nlohmann::json geojson;
    
    geojson["type"] = "FeatureCollection";

    geojson["name"] = "Eddy_Analysis_" + method_name + "_" + std::to_string(static_cast<int>(depth_m)) + "m";
    
    // Set coordinate system
    nlohmann::json crs;
    crs["type"] = "name";
    nlohmann::json crs_properties;
    crs_properties["name"] = "urn:ogc:def:crs:OGC:1.3:CRS84";
    crs["properties"] = crs_properties;
    geojson["crs"] = crs;
    
    // Create a feature array
    nlohmann::json features = nlohmann::json::array();
    
    for (const auto& eddy : eddies) {
        features.push_back(createFeature(eddy, method_name, depth_m));
    }
    
    geojson["features"] = features;
    
    return geojson;
}

std::string GeoJSONExporter::exportToGeoJSON(const std::vector<EddyContour>& eddies,
                                              const std::string& method_name,
                                              double depth_m) {
    nlohmann::json geojson = createGeoJSONObject(eddies, method_name, depth_m);
    
    // Convert to a formatted JSON string
    return geojson.dump(2);
}

nlohmann::json GeoJSONExporter::create3DEddyFeature(const Eddy3D& eddy_3d) {
    nlohmann::json feature;
    feature["type"] = "Feature";

    // Geometry uses the boundary of surface vortices
    if (eddy_3d.layers.empty()) {
        feature["geometry"] = nullptr;
    } else {
        const auto& surface_layer = eddy_3d.layers.front();
        nlohmann::json ring = nlohmann::json::array();
        for (const auto& point : surface_layer.boundary_points) {
            ring.push_back({point[0], point[1]});
        }
        // Ensure closure
        if (!ring.empty() && ring.front() != ring.back()) {
            ring.push_back(ring.front());
        }
        feature["geometry"] = {
            {"type", "Polygon"},
            {"coordinates", {ring}}
        };
    }
    
    // Rich 3D attributes
    nlohmann::json properties;
    properties["id"] = eddy_3d.id;
    properties["type"] = eddy_3d.type;
    properties["layer_count"] = eddy_3d.layer_count;
    properties["volume_km3"] = eddy_3d.volume_km3;
    properties["center_lat"] = eddy_3d.center_3d[0];
    properties["center_lon"] = eddy_3d.center_3d[1];
    properties["mean_depth"] = eddy_3d.center_3d[2];
    properties["min_depth"] = eddy_3d.vertical_extent[0];
    properties["max_depth"] = eddy_3d.vertical_extent[1];
    
    // Serialize the information of all layers
    nlohmann::json layers_array = nlohmann::json::array();
    for (const auto& layer : eddy_3d.layers) {
        layers_array.push_back({
            {"depth", layer.depth},
            {"area_km2", layer.area_km2},
            {"mean_vorticity", layer.mean_vorticity},
            {"mean_omega", layer.mean_omega},
            {"center_lat", layer.center[0]},
            {"center_lon", layer.center[1]}
        });
    }
    properties["layers"] = layers_array;

    feature["properties"] = properties;
    return feature;
}

std::string GeoJSONExporter::export3DEddiesToGeoJSON(const std::vector<Eddy3D>& eddies_3d,
                                                       double depth_m) {
    nlohmann::json geojson;
    geojson["type"] = "FeatureCollection";
    
    nlohmann::json features = nlohmann::json::array();
    for (const auto& eddy : eddies_3d) {
        features.push_back(create3DEddyFeature(eddy));
    }
    geojson["features"] = features;
    
    return geojson.dump(2);
}

} // namespace EddyDetection 